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^ Abstract 

d We introduce an interpolation framework for T^oo niodel reduction founded on ideas 

H originating in optimal-7^2 interpolatory model reduction, realization theory, and com- 

'~~' plex Chebyshev approximation. By employing a Loewner "data-driven" framework 

^vq within each optimization cycle, large-scale T-Loo norm calculations can be completely 

!>■ avoided. Thus, we are able to formulate a method that remains effective in large-scale 

settings with the main cost dominated by sparse linear solves. Several numerical exam- 

(^ pies illustrate that our approach will produce high fidelity reduced models consistently 

^O exhibiting better Tioo performance than those produced by balanced truncation; these 

C^ models often are as good as (and occasionally better than) those models produced by 

optimal Hankel norm approximation. In all cases, these reduced models are produced 

at far lower cost than is possible either with balanced truncation or optimal Hankel 

norm approximation. 



X 

"^ 1 Introduction 

The need for high accuracy mathematical models in problems involving simulation and con- 
trol often results in dynamical systems described by a large number of differential equations. 
Working with such large-scale systems can easily place overwhelming demands on compu- 
tational resources, a problem which model reduction methods seek to alleviate by approxi- 
mating the original model with another model consisting of far fewer (but carefully crafted) 
differential equations. Strategies for carrying out this approximation should be both efficient 
and accurate. For an overview of model reduction, see [Ij. 



We consider here single-input/single-output (SISO) linear dynamical systems given in state- 
space form as: 

Ex{t) = Ax{t) + bu{t), y{t) = c^x{t) + du{t), (1) 

where E, A & M"^", b,c & M" and d E R. E is assumed to be nonsingular throughout, 
although our approach extends without difficulty to cases where E is singular so long as 
nullity(^) = nullity(^^) (that is, provided that is a nondefective eigenvalue for E). x{t) G 
M" are the states] u{t) G M is the input] and y{t) G M is the output of the dynamical system 
in ([l]). The transfer function of the system is 

H{s) = c^ {sE - A)-^b + d, 

defined for s G C. In accord with standard convention, we denote both the system and 
its transfer function by H{s). We assume that H{s) is both controllable and observable. 
The order of H{s) is the number of poles it possesses, counting multiplicity. Since E is 
nonsingular, all poles of H{s) are finite and since H{s) is both controllable and observable, 
the order of H{s) is identical to the dimension, n, of the state vector a; in ([I]). 

We denote by "H^, the set of rational functions of order at most k which are bounded and 
analytic in the closed right half plane in C. We assume in all that follows that H G "HJ^. 
The "Hoo norm of H is defined as 

Ili^IlT, = max I Hiioj) I . (2) 

1 

2 
|2 



If the input function, u(t), is square integrable: ||w||l2 = I / | M(i) | dt\ < oo, then 

the output function, y(t), of ([I]) will be square integrable as well; u and y have Fourier 
transforms u, y E L^(m) that are related according to y{u)) = H{juj)u{uj). One immediately 
observes that the "Hoo norm defined in d2| is an L^-induced operator norm of the underlying 
convolution operator mapping m H- y. 

Our goal is to construct another system 

ErXr{t) = ArXr{t) + bj.u{t), yr{t) = cjxr{t) + drU{t) (3) 

of much smaller order r -^ n, with Er, A^ G M''^'", bj., c^ G M*", and d.,. G M determined so 
that yr approximates y uniformly well over all u G L^(E"'"), in an appropriate sense. 

Toward this end, define a reduced transfer function associated with pi) as Hr{s) = cJ(sEr — 
Ar^^br + dr. Then, for any u G L2(M+), 

Wy-yrh^ < \\H - Hr\\noo\\u\\L2. (4) 

The error transfer function, H{s) — Hr{s), may be rewritten as 

H{s) - Hr{s) = c^{sE - A)-^b - [cJ{sEr - Ar)-% + {dr - d)] . 



Evidently, a non-zero (i-term in the original model may be absorbed into a reduced-order 
model by assigning dj. H- (rf^ — d). This allows us to assume in all that follows that d = in 
the original model without any loss of generality. 

In view of Q, the output error magnitude, \\y — yrWi'^, may be made uniformly small over 
all bounded inputs u (say, with ||m||l2 < 1) if we find a reduced system, Hr, that makes 
\\H — HrW-Hoc small. This leads naturally to the optimal Hoo model reduction problem: 



For a given H G Ti'^ and reduction order r < n, find H* G T-L^, that solves 



mm 



H — Hr 



rio 



(5) 



This problem is an active area of research [2] . We develop a methodology for approximating 
solutions to (|5| that remain effective in large-scale settings, settings where the original state- 
space dimension, n, could be on the order of 10^ or more, for example. Most methods known 
to us will be intractable even for modest system order, say, on the order of a few thousand 
(with the exception of balanced truncation, see the discussion below). 

Kavranoglu and Bettayeb [26j showed that (tsl) can be converted into an optimal Hankel 
norm approximation problem for a special imbedded system with augmented input and 
output mappings. However, this approach is infeasible in practice since knowledge both 
of the minimum of (tsl) as well the imbedded system is required. As noted in [21], this 
information is available (or computationally accessible) only in very special cases. 

Several methods to solve (IS]) that utilize linear matrix inequality (LMI) frameworks have 
been presented as well; see, for example, ^5\ |25l |2l[ [231 E] and references therein. These 
approaches rapidly become computationally intractable with increasing state space dimen- 
sion. Indeed, published examples illustrating LMI-based methods in [TSl [2S1 [211 [231 HI] all 
had order less than n = 10. 

The most common practical methods for obtaining satisfactory "Hoo reduced models are 
Gramian-based methods such as balanced truncation (BT) [311 [32] and optimal Hankel norm 
approximation methods (HNA) [H]. Both approaches are known to yield small approximation 
errors in the "Hoo norm [121 H] , though neither generally is capable of producing globally opti- 
mal solutions to (|5|. Both approaches also remain computationally feasible for modest state 
space dimension (perhaps a few thousand), but significantly larger state space dimension 
still presents challenges. HNA requires an all-pass dilation of the full-order model followed 
by a full eigenvalue decomposition. These are dense matrix operations, having a complexity 
growing with 0{n^), which generally limits problem sizes to a few thousand. Notably, [TT] 
was able to extend HNA to dynamical systems having state space dimension of 0(10^), using 
state-of-the-art numerical techniques tailored to high performance computer architectures. 

The situation for BT is a bit better. BT has been applied to systems of order (9(10^) by 
solving the underlying Lyapunov equations iteratively using ADI-type algorithms; see, for 
example, [2Q1 [381 ESI [lOl [351 122] and references therein. 



We describe here a different model reduction methodology that can be applied effectively 
even for very large state space dimension, yielding reduced models typically having smaller 
T-Loo errors than either BT or HNA moreover, at significantly lower cost. Towards this end, 
we present a new framework for the "Hoc approximation problem using interpolatory model 
reduction. By connecting ideas from interpolatory I-L2 model reduction [21], realization 
theory [29], and complex Chebyshev approximation [39], we develop an interpolation-based 
method for 'Hoo-approximation that remains numerically efficient in large-scale settings. The 
main cost of our approach involves the solution of sparse linear systems. We demonstrate 
that our approach can yield reduced-order systems having l-L^ errors which are often half 
that of BT, and very close to (indeed, sometimes better than) that of HNA. For symmetric 
systems, our method typically produces reduced-order systems with "Hoo-errors that are near 
the theoretical best possible of (|5]). 

The rest of the paper is organized as follows: we close this section with a brief review of 
interpolatory model reduction and related approaches for solving the optimal 7^2 approxi- 
mation problem. We introduce our method in ^ and illustrate its effectiveness via several 
numerical examples in ^ 



Interpolatory model reduction. Given a dynamical system H[s) and a set of points 
{si, S2, . . . , Sk] C C, interpolatory model reduction produces a dynamical system Hr{s) 
such that Hj.{s) interpolates H{s) together with a prescribed number of derivatives at the 
points {si, S2, ..., Sk}. Although this is posed as a rational interpolation problem, the 
construction of a solution may be accomplished with a variety of rational Krylov subspace 
projection techniques. Rational interpolation via projection was first proposed by Skelton 
et al. [m HH US] . Later, Grimme [16] showed how to construct a reduced-order (Hermite) 
interpolant using a method of Ruhe. 

Theorem 1.1 (Grimme [IE]). Given H{s) = c^{sE — A)^^b and a point-set Si G C 
containing r distinct points: Si = {si, . . . , Sr}, let 



Vr = [{siE - A)-^b . . . {srE - A)-^b] 



Wj 



'c^{siE - A)-^' 



c^isrE - A)-^ 



(6) 



Define a reduced-order model H^.{^) — c^isEj. — Aj.) ^br, where dr = and 
Er = WjEVr, Ar = WjAVr, 6, = W^b, and c^ = c^Vr. 



(7) 



Then H{si) = H'^isi) and H'{sk) = H^'{sk), for i = 1, . . .2r where ' denotes the derivative 
with respect to the frequency parameter, s. 



Higher-order derivatives can be matched similarly; for details, see [T^ 1^. 



'H2-optimality conditions. Theorem 1A_ gives explicitly computable conditions that will 
yield a reduced-order model satisfying 2r interpolation conditions; one need only solve 2r 
linear algebraic systems to form the columns of Vr and Wr- Notably, Theorem 1.1 carries no 
hint of how best to choose these interpolation points. A method for determining interpolation 
points that leads to reduced-order models that are (locally) optimal with respect to the H2 
error was developed in |21j and will be a point of departure for the approach we propose 
here. 

For a SISO dynamical system H{s), the 7^2 norm is defined as 



1 



H\\n,= [i^ I I HM \' du 



1/2 



Then, for a given full-order model H{s), and selected reduction order r < n, the optimal ^2 
approximation problem seeks a reduced model, H^{s), that solves 



mm 



H-m 



(9) 
■H2 



The approximation problem ^ has been studied extensively; see, for example, [30] . |13] . [21] . 
[36] . [iO] . [T7] . [6], [Z],[1Z] and references therein. First-order necessary conditions for "^2- 
optimal approximation may be formulated in terms of interpolation conditions: 

Theorem 1.2 ([30l[2T]). Given a full- order model H{s), let H^{s) be an 'H2-optimal reduced 
order model of order r, with simple poles \i, . . . , X^. Then, 

i7(-A,) = /f°(-Ai) and H'{-X,) = H^'{-X,) for 2 = l,...,r, (10) 

where ' denotes differentiation with respect to the frequency parameter, s. 

These necessary conditions characterize the 'H2-optimal reduced order model as a rational 
Hermite interpolant matching the full-order transfer function and its derivative at mirror 
images of the reduced-order system poles. This can be accomplished with the help of Theo- 
once the poles of H^{s) are known. Of course, the poles of -ff°(s) are not known a 



1.1 



rem 

priori, but they can be computed iteratively using the Iterative Rational Krylov Algorithm 

(IRKA) developed by Gugercin et al. [21] . 

IRKA is a fixed point iteration that in the SISO case typically exhibits rapid convergence to 
a local minimizer of the ?^2-optimal model reduction problem. Sparsity in E and A can be 
well-exploited in the linear solves of Steps 2 and 3c and IRKA has been remarkably successful 
in producing high fidelity reduced-order approximations in large-scale settings; it has been 
applied successfully in finding 'H2-optimal reduced models for systems of high order (e.g., 
n > 160, 000, see [27]). For details on the algorithm, we refer to the original source [2T] . 



Algorithm IRK A. Iterative Rational Krylov Algorithm J21^ 

Given a full-order H{s), a reduction order r, and convergence tolerance tol. 

1. Make an initial selection of interpolation points Sj, for i = 1, . . . , r that is 
closed under complex conjugation. 

2. Construct Vr and W^ as in ^ with Sj+r = Sj for z = 1, . . . , r 

3. while (relative change in {sj} < tol) 

a) Er = W^EVr and A, = W^AV, 

b) Solve r X r eigenvalue problem A^u = A-E^u. 

c) Assign Si = Sj+r < Xi{Ar, Er) for i = 1, . . . , r. 

d) Update V,. and Wr as in ^ using new {si}. 

4. Er = WjEVr, Ar = W^AVr, K = W^b, c^ = c^Vr, and dr = 0. 



2 An interpolatory approach for Hoo approximation 

The principal result that defines the character of our approach was provided by Trefethen 
in [3n]; it is an analog to the Chebyshev Equioscillation Theorem. 

Theorem 2.1 (Trefethen [39]). Suppose H{s) is a transfer function associated with a dynam- 
ical system as in pi). Let m^^^^^s) be an optimal "Hoo approximation to H{s) (i.e., a solution 
to ^ ) and let Hr{s) be any r^^ order stable approximation to H{s) that interpolates H{s) 
at 2r -\- 1 points in the open right half plane. Then 



min \H{ju) - Hr{juj)\ < \\H - i^riko. < 11^ " Hr 



Wa 



ojeJi 



In particular, if \H{]uj) — Hr{]U))\ = const for all u E R then Hr{s) is itself an optimal "Hoo 
approximation to II{s). 

One sees from this that a good Tioo approximation will be obtained when the modulus 
of the error, \H{s) — Hr{s)\, is nearly constant as s = ju runs along the imaginary axis. 
We select 2r + 1 interpolation points in the open right half plane that will induce this. By 



utilizing Theorem 1.1 , we may locate 2r interpolation points in the right half plane as we like, 
producing an interpolating reduced-order system, -ff°(s). Also, H{oo) = if°(oo) = so we 
can exploit the freedom in choosing dr to move the (2r + 1)*^* interpolation point from oo into 
the open right half-plane. Note that the straightforward construction, Hr{s, dr) = H^{s)+dr, 
creates a reduced-order model that has all interpolation points depending on dr- We prefer a 
different formulation that uncouples the 2r interpolation points associated with H^{s) from 

6 



the influence of dr- Tlie construction that acconiphshes this was introduced in [29] (see also 
[5] for a formulation close to what we use here). 

Theorem 2.2. Given H{s) = c^{sl — A)^^h and a point-set S G C containing r distinct 
points: S = {si, . . . , Sr}, let Vr, Wr, Ej., Ay, br, and Cr be defined as in Theorem \l.l\ For 
any given dr G M, define a new reduced- order system 

Hr{s, dr) = {Cr — drSf {sEr — Ar — dree^)^^{br — drC) + dr (11) 

where e denotes a vector of ones. Define auxiliary reduced systems: 

iyO(s) = c^{sEr - Ary^br, Gi{s) = e^{sEr " A,)-i^, 
G2{s) = cJ{sEr - Ar^^e, and (^3(5) = e^{sEr - ArY^e. 

Then 



^0.„^ , . (Gi(s)-l)(G2(s)-l) 



Hr{s, dr) = H'r{s) + dr ' '^ [ - dG^) ^^^^ 

and for all dr G M 

H{s,) = H^{si)=Hr{si,dr) and H'{s,) = H^/{s,) = Hl{s,,dr) for t = l,...,r 
where ' denotes the derivative with respect to the frequency parameter, s. 



Proof: The expression (12) follows from (11) with straightforward manipulations that begin 



with the Sherman-Morrison formula: 

{sEr - Ar - dree^)-^ = {sEr - Ar)'^ + /' {sEr - Ar)'^ee'^ {sEr - Ar)'^ 

1 — drGsls) 

Define P{s) = Vr{sEr — Ar)^^W^ {sE — A) . Observe that P{s) is a (skew) projection onto 
Ran(l^) for any s G C for which it is well-defined and so we have 

VrBk =PiSk)Vrek = [VriSkEr - Ar)-'W^iSkE - A)] {SkE - A)-^b 
= Vr{SkEr - Ar)'^W^b = Vr{SkEr - Ar)~^br 

Linear independence of the columns of Vr then implies e^ = {skEr — Ar)~^br and thus, 
Gi{sk) = 1, for A; = 1, 2, . . . , r. A similar argument yields G2{sk) = 1, for k = 1, 2, . . . , r. 
Taken together, we get that H'^{si) = Hr{si, dr) for i = 1, 2, . . . , r. Likewise, 

H'(.ri) ffO'r,^ - ^^^ki^Lrr r,^ run.) 1^ , , g;(^)(g2(.)-i) , , jCis) - i)G'2is) 

HAs,dr)-Hr is) - ^^ _ d^G^^,^)2(Gl{s)-l){G2{s)-l)+dr ^ _ ^^^^(^) +^r ^ _ ^^,^,{3) 

SO that H^'{si) = H'^{si,dr) as well. D 

Let {Hr{s, ■)}dr denote the set of all transfer functions Hr{s, dr) with dr ranging over M. The 
freedom we have in choosing dr is significant to us for at two reasons. First, {Hr{s, ■)}d^ is a 

7 



parameterization of the set of all proper rational functions of degree r having real coefficients 
that satisfy the same interpolation constraints as -ff°(s) (see e.g., [29] )• Second, it is now 
possible to construct reduced-order models of order r satsfying 2r+l interpolation conditions, 
which is an essential step towards constructing reduced-order models that are optimal in the 
Hoo-norm. Since Hr{s,dr) interpolates H{s) at si,...,S2r for any dr, one could select an 



additional (real) interpolation point, S2r+i > and directly calculate from (12) the value of 
dr that enforces Hr{s2r+i,dr) = H{s2r+i)'- 

, ^ H{s2r+l) - H^{s2r+l) 

' (Gi(s2.+ l) - 1)(G2(S2.+1) - 1) + Gsis2r+l)iHis2r+l) - i/°(s2.+l)) ' 

We avoid the necessity of explicitly selecting S2r+i- Instead, as discussed below, dr will be 
chosen directly to decrease the 7/oo error. 



2.1 An algorithm for 7/oo approximation 

It has been observed (e.g., see [2T],|1]) that 7^2 optimal interpolation points produced by 
IRKA yield reduced models that are not only (locally) 'H2-optimal but frequently also are 
high-fidelity Hoc approximations to the original system. Indeed, H2-optimal models pro- 
duced by IRKA yield 7/oo error norms that are comparable to that of BT and sometimes are 
even better. Therefore, our approach begins with Algorithm IRKA to obtain 2r interpolation 
points (counting multiplicity - Hermite interpolation at r distinct points) determining an 
'H2-optimal reduced model, H^{s). This choice for H^{s) defines a family of approximations 
parameterized hy dr, {Hr{s,-)}dr- We then proceed by (approximately) minimizing the "Hoo 
error with respect variations in dr. These steps are summarized below: 



Algorithm I HA. Interpolatory "Hoo Approximation: 

Given a full-order model, H{s), and reduction order, r. 

1. Apply Algorithm IRKA to compute 2r ?^2-optimal interpolation points and an 
associated 'H2-optimal reduced model, H^{s). 



2. Find dt = argmin \\H — i/r. IL where Hr = Hr(s,dr) is defined in (12) 



3. Construct the final Hoo approximant as H*{s) = Hr{s,d*). 



The distribution of interpolation points obtained in Step 1 (as an outcome of 'H2-optimal 
approximation) yields very effective Hoo approximants as well. We therefore wish to preserve 



these interpolation points using Theorem 2.2 while varing the dr parameter in such a way 



as to drive down the "Hoo error - centering the error curve about the origin in the process. 
H*{s) will denote an "Hoo approximant having 1) an 7^2 optimal pole distribution (from Step 
1 of Algorithm IHA) and 2) an optimally chosen d* (from Step 2). 



2.2 Efficient Implementation of Step 2 of Algorithm IHA 



The major contributions to the cost of IHA come from hnear solves arising in Step 1 (from 
IRKA) and the "Hoo norm evaluations required in solving the (scalar) nonlinear optimization 
problem in Step 2. Hoo norm evaluation involves repeated solution of several large-scale 
Riccati equations of order n + r. Solving even a single Riccati equation, let alone several, 
will be a formidable task when n is on the scale of tens of thousands or larger, the range 
of system dimension of interest here. We describe below an effective strategy to circumvent 
this difficulty. 



The optimization problem of Step 2 can be rewritten (from Theorem 2.2) as 



mm 

drGR 



H{s) 



H^( 



dr{G,{s) - 1){G2{S) 



1 - drG^is) 



flo 



where H^{s) is a reduced-model obtained by IRKA in Step 1 of Algorithm IHA. 

If one can find a reduced-order approximation, Fk{s), to the error system, F[s) = H{s) — 
H^{s), having modest fidelity and order k <t^ n, then an associated optimal dr-term. could 
be efficiently calculated by solving the (comparatively) low order optimization problem 



mm 

d,.eiR 



F,is) 



dr{G,{s)-l){G2{s)-l) 



1 - drG^is) 



(13) 



Ha 



Provided k <t^ n, the cost of solving (13) will be negligible compared to that of original 
problem. Of course, whatever advantage this strategy may bring could be nullified if the 
cost of obtaining Fk{s) is significant. By using a Loewner matrix approach developed by 
Mayo and Antoulas p9] and described briefiy below, we are able to obtain Fk{s) at negligible 
cost relative to the computational demands already incurred in Step 1. We reuse information 
obtained in the course of IRKA in Step 1 to obtain, for negligible additional effort, a reduced 
error model Fk{s) having modest fidelity, adequate for the demands of Step 2. 

Suppose we have evaluated the error system, F{s), and derivative, F'{s), on a set of distinct 
points {si, S2, . . . , Si} C C. We will construct from this data a reduced order surrogate. 



Fk(s) 



[sEk - Ai 



^bh so that 



F{si) = Fkisi) and F'{si) = F^(si) for ^ = 1,2, 



The Loewner matrix approach as developed by Mayo and Antoulas [29] permits "data- 
driven" model reduction; one need not have access to state-space matrices determining a 
realization of the full order system. Only "response measurements" are used, that is, transfer 
function evaluations. Reduced-order models will be constructed directly that interpolate this 
"measured data". 



Define matrices L G C' 


^^ and M e d^^ 


as 


(L). . := 


r F{s^) - 

< 


^(^.) 
^i 


if Zy^J 

if z = j 





SiF{si) - SjF{sj] 



b{ bj 



if i^j 

^M- , (14) 

[.F(.)]'|_^ if ^=j 

L is the Loewner matrix associated with interpolation points si,S2, ■ ■ ■ ,se and the dynamical 
system -F(s); M is the corresponding shifted Loewner matrix (see [22] for details). Once L 
and M are constructed, assume that the interpolation data satisfy the following assumption: 

L 



rank (s^L - M) = rank[L M] = rank 



(15) 



for i = 1,2, . . . ,i. For SISO systems, this assumption holds whenever an interpolant of order 
r = rank (sjL — M) exists [29]. For MIMO systems, this assumption is also generically valid, 
but details are more involved; the interested reader should see Lemma 5.4 of 129). 



rank 



In light of (|15|), a rational Hermite interpolant is constructed by first choosing k so that 

L 



> k. Then for some choice oi 1 < i < i, compute SjL — M = Y@X*, the SVD 



of SjL — M. Let Yfc G c^^^ and X^ G c^^'^ denote the leading k columns of Y and X, respec- 
tively (associated with a truncated SVD of order k). Let Z = [F{si), F{s2) , . • • , F{si)]'^ 
and define 

Ek = -Y:hXk, Ak = -Y:MXk, bk = Y:Z, Ck = Z^Xk, 

and Ffc(s) = cj^i^sEk — Ak)~^bk. k may be considered as a truncation index here. Depending 
on whether k is chosen so that k = rank (sjL — M) or k < rank (sjL — M), Fk{s) will then be 
either an exact or an approximate interpolant, respectively. In practice, one should choose k 
no larger than the numerical rank of SjL — M, which can be determined by a Singular Value 
Decomposition (SVD) of SjL — M. See [21] for a full development of these ideas. 

Observe that until convergence occurs within Step 1 of Algorithm I HA, every cycle of Step 3 in 
Algorithm IRKA will provide a sampling of H{s) and H'{s) at r interpolation points - generally 
a different set of interpolation points in each cycle. (Step 3 of Algorithm IRKA constructs a 
Hermite interpolant on a set of r interpolation points that is cyclically adjusted.) Suppose 
that Algorithm IRKA takes q steps to convergence. When Step 1 of Algorithm IHA concludes, 
we will have had H{s) and H'{s) sampled at a total oi i = q x r interpolation points. We 
collect these interpolation points and transfer function evaluations throughout IRKA. Once 
IRKA converges, (yielding an 7^2-optimal model, H^{s)), we evaluate if° and H^' at these 
i points as well. Since the order of H^{s) is r, the cost of these function evaluations is 
negligible. After the completion of Step 1 of Algorithm IHA, we have an £-iold sampling 
of both F{s) = H{s) — H^{s) and F'{s) with virtually no additional computational cost 
beyond what was needed for Step 1 itself. Then, we simply apply the Loewner matrix 
approach described above to construct Fk{s). The choice of k will be clarified via numerical 
examples in ^ 
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Numerical Cost of IHA. Note that once Step 1 of Algorithm IHA is completed, Step 2 is 
not computationally intense. The cost is dominated typically by an £ x £ SVD computation 
with i = q X r. Since IRKA typically converges rather quickly (especially so in the SISO 
case focused on here), i is generally modest in size. In all of our numerical examples, we 
have never needed to compute an SVD of size larger than 200 x 200; a trivial computation. 
Moreover, in all of our numerical examples k never exceeded 33; making the solution of the 
optimization problem in Step 2 quite cheap. The overall cost of IHA is only marginally more 
than that of IRKA and is dominated by the same sequence of sparse linear solves. 



Stability of the reduced model: The asymptotic stability of Hj. in Step 2 may be 
enforced by adding a penalty function to the cost function penalizing values of d^ that yield 
systems having poles too close to the imaginary axis. The optimization algorithm would then 
automatically reject dr terms that cause unstable eigenvalues in the pencil sEr — Ar — drCe^ . 
Alternatively, one may simply calculate and check the eigenvalues of this r x r pencil to 
determine whether to accept a dr on the basis of stability. In our numerical examples, we 
simply reject dr values that cause unstable reduced models by setting the corresponding the 
function value to oo. We always obtained a stable reduced- model as a result but there are 
better, more effective numerical strategies to perform this task. For example, a logarithmic 
barrier function that takes the real part of the pole closest to the imaginary axis as its 
argument. These numerical issues will be studied in a separate work where we extend the 
method to the MIMO as discussed next. 



Application to MIMO systems: Algorithm IHA, together with the results of Section 2.2, 
can be easily generalized to the MIMO case. In the MIMO case, IRKA enforces bitangential 
interpolation conditions at the reflection of the reduced-order poles across the imaginary 
axis. A complete account of 7^2 optimal model reduction and IRKA for MIMO systems can 
be found in [4j. These interpolation conditions can be enforced while varying the D-term 
of the reduced-order system. See Theorem 3 in [5] for a complete description of how to 
construct the bitangential interpolant while varying the D-term; thus the theory in this 
paper directly generalizes to the MIMO case. The optimization step involving the matrix 
D G IRP^™-, however, is more involved than its scalar counterpart for SISO systems and a 
robust numerical implementation is the subject of ongoing research. 



3 Numerical Experiments 

We illustrate here the performance of IHA on various benchmark models for model reduction 
and compare its performance with that of balanced truncation (BT) optimal Hankel norm 
approximation (HNA), and Iterative Rational Krylov Algorithm (IRKA). We note that generic 
BT yields dr = 0, a. null feed-forward term. Therefore, in order to present a fair representation 
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for BT, we add a final step that varies the dr term in the BT model as well; we use d^ that gives 
the minimum "Hoo norm. We refer to BT with this optimally chosen d^ term as "modified 
balanced truncation" for MBT). 



3.1 PEEC Model 

The full-order system is the spiral inductor system PEEC model \^2\ of order n = 1434. The 
system is state-space-symmetric (SSS), i.e. the transfer function H{s) = c^{sE — A)~^b 
satisfies E = E^ > 0, A = A^ and c = b. SSS systems appear in many important 
applications such as in the analysis of RC circuits, and has been the subject of several 
model reduction papers; e.g, see [28l [371 [3ll HB]. We first illustrate the effect of the dr-teim. 
modification for MBT and for Step 2 of Algorithm I HA. In MBT, once the initial BT phase is 
completed, we vary dr and measure the resulting changes to the "Hoc error. For I HA, once the 
IRKA phase in Step 1 is completed, we vary the dr term (inducing corresponding changes to 



Ar, hr and c^ as in Theorem 2.2) and once again measure the resulting changes to the T-Loo 
error. Note that calculation of T^oo error is for illustration only and is not used by IHA to 
compute the optimal dr-term. Results are shown in Figure [T] for r = 2. As the figure shows, 
there is essentially no improvement in MBT that comes from adjusting dr, however for IHA 
the error is reduced by a factor larger than two. Even though at the starting point, dr = 0, 
the 'H2-optimal approximation has higher "Hoc error than does BT at dr = 0, IHA is able 
to reduce the "Hoc error to a value significantly lower than that for MBT through dr-teim 
optimization, at virtually negligible computational cost. In Figure [T| the point dr = on the 
curve for IHA gives the value of the "Hoc error produced by IRKA. Note that the "Hoo error for 
IHA is less than half of that for IRKA. This behavior is common to all the numerical examples 
that follow. 

Before presenting the comparisons between IHA, MBT, and HNA, we illustrate the efficiency 
of the methodology outlined in §2.2 in solving the optimization problem in Algorithm IHA; in 



other words in finding the optimal dr-term in Figure [TJ For the three r values r = 2, 4, 6, we 
implement IHA both by exactly solving Step 2 and by the method of ^2.2 Table [T] tabulates 



the results where the resulting optimal dr values and the T^oo error norms for both methods 



together with the order-fc used in method of ^2.2 are listed. As the table clearly illustrates 



that the Loewner matrix approach yields dr terms and the "Hoo error norms which are very 
close to true-optimal values of the underlying optimization problem in Step 2 of Algorithm IHA. 
More importantly, this is achieved with negligible computational cost where the function 
evaluations are the "Hoo norm computations for an order k system only as opposed to order 
n + r. One pattern we have observed via several SSS models is that A; = 2r -|- 1 is a natural 
choice. This pattern repeated itself for every SSS example we have tried. There was a clear 
cut-off point in the singular values of the matrix SjL — M at (2r + l)*"^ singular value. The 
decay of these singular for all three r values are shown in Figure [2] supporting the k = 2r + 1 
choice. 
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Figure 1: Comparison of the "Hoo Error as the dr-term varies for IHA and MBT 



Table 1: Solution of the optimization problem in Step 2 
Exact Loewner 



dl 



\H — Ht 



l-Ha 



dt 



\H — H* 



I Wo 



6.9577 X 10-3 
1.0041 X 10"^ 
2.7795 X 10-6 



4.4522 X 10-3 
8.6577 X 10-5 
4.4771 X 10-6 



6.9659 X 10-3 4.4574 x lO'^ 5 
1.0076 X lO-'^ 8.7114 x 10-^ 9 
2.7804 X 10-6 4.4857 x 10-^ 13 



Next, to compare IHA, MBT, and HNA, we reduce the order of the system to r = 2,4,6. 
The resulting relative Tioo error values together with the lower bound (i.e. o"r+i/ |-f^||-Hoo 
where cxr+i denotes the (r + 1)**^ Hankel singular value of H(s)) are listed in Table 2) The 
lowest error value for each r is shown in bold font. For every r, IHA outperforms MBT by 
almost a factor of two. Also it performs very close to HNA, indeed outperforms it for r = 4. 
This shows the strength of the proposed method. Without solving any Lyapunov equations 
and without the need for any large-scale "Hoc norm computation, our methods consistently 
outperforms BT by a significant amount and performs as nearly as and sometimes better 
than HNA. Note that T-Loo error values for IHA is very close to the lower bound given by 
ar+i/WHW-n^. We can relate this to Theorem 2.1 for each order of approximation shown in 
Table [2| the reduced-model due to IHA results in exactly 2r + 1 interpolation points in the 
right-half plane and a nearly circular error curve as illustrated in Figure [3] . 2r of these zeros 
result from IRKA; indeed these are r distinct zeroes with multiplicity 2. Then, the (2r -|- 1)**^ 
zero is obtained by introducing the dr term. For example, for the case r = 6, computing the 
optimal dr-term in turn placed an additional zero at the point 6.60 x 10*^. 
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Normalized singular values of s L - M 




10 12 14 16 18 20 



Figure 2: The decay of the singular values of SjL — M for PEEC Model 
Table 2: Relative "Hoo error norms for the PEEC Model 



r 


IHA 


MBT 


HNA 


Lower bound 


2 


4.45 X 10-=^ 


8.21 X 10-3 


3.95 X 10 =^ 


3.72 X 10-3 


4 


8.66 X 10^ 


2.78 X 10-^ 


1.24 X 10-^ 


7.79 X 10-5 


6 


4.48 X 10-6 


8.16 X 10-6 


3.41 X 10^ 


3.15 X 10-6 



Note that the findings in this SSS example are common to all other SSS examples we have 
tried. The T-Lao error due to IHA has always been around half of MBT, and fc = 2r + 1 has 
been a clear cut-off point for the Loewner matrix approach. Due to space limitations, we 
omit these examples. For some other SSS examples, we refer the reader to J13j . 



3.2 CD Player Model 

Next, we demonstrate the proposed method on the CD player model of order n = 120. For 
details on this model, see [31 [12]. We reduce the order to r = 2, 4, 6, 8, 10 using IHA; stopping 
at r = 10 as the relative error fell below 10-^. As for the previous example, we compare 
the effect of solving the optimization problem in Step 2 of Algorithm IHA exactly versus by 
the method of §2.2[ The results are listed in Table [3] As for the previous case, the Loewner 
matrix approach yields 'Hoo error norms which are very close to the true optimal values of 
the underlying optimization problem. Once again, the function evaluations are much simpler 
as k never exceeds 33. In this example, we chose the value of k as the normalized singular 
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Figure 3: The Nyquist plot of the error system H — H* for PEEC Model 



values of SjL — M drops below the tolerance of 10 ^. Figure |4] shows this decay behavior for 
r = 4 and r = 6. 



Table 3: Solution of the optimization problem in Step 2 
Exact Loewner 



dt 



\H — Ht 



\V.c 



\H - m 



I lio 



2 
4 
6 
8 
10 



-3.6728 
3.5019 X 10-1 
2.9538 X 10-1 
1.3888 X 10-1 

-3.5750 X 10-2 



3.6597 X 10-1 
2.1318 X 10-2 
1.0155 X 10-2 
4.8357 X 10-3 
8.5384 X 10-^ 



-3.6545 
2.4819 X 10-1 
2.0763 X 10-1 
1.3625 X 10-1 

-3.2438 X 10-2 



3.6604 X 10-1 2 

2.1422 X 10-2 20 

1.0426 X 10-2 25 

4.8526 X 10-3 32 



8.9952 X 10 



-4 



33 



Next, we compare IHA with MBT and HNA. The results are illustrated in Table |4] where the 
minumum value for each r is shown in bold font. Note that for every r value, the proposed 
approach outperforms MBT. Moreover, except for the r = 2 and r = 4 cases, IHA outperforms 
HNA as well. 

Before moving to the next example, we illustrate the effect of the (i^-term modification in 
our proposed method as opposed to MBT. In Figure [5} we show how the absolute "Hoo-error 
changes both in IHA and in MBT as we vary the dr-term for r = 10. In both cases, (ir = is the 
starting point. While the ?^oo-error reduces marginally from 0.0905 to 0.0852 in MBT- only 
a 5.86% reduction, the gain is much more significant in IHA where we reduce the "Hoo-error 
from 0.0938 down to 0.0585, a significant 37.67% reduction in the Hoo error. Even though 
the 'H2-optimal approximation has a larger l-Loo error at the starting point, (i^ = 0, than does 
the BT approximation, the (ij,-term optimization in IHA will produce a final "Hoo error that is 
significantly lower than that produced by the corresponding (i^-term optimization in MBT. 
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Normalized singular values of s L - M 





Figure 4: The decay of the singular values of SiL — M for the CD Player Model 

Table 4: Relative "Hoo-norm Error norms for the CD Player Model 

r I IHA MBT HNA 

^ 3.66 X 10-^ 3.68 x lO^^ 3.35 x lO^^ 
4 2.14 X 10-2 2.25 x 10"^ 2.00 x lO'^ 
6 1.04 x 10-2 1.19 X 10-2 1.23 x 10"^ 
8 4.85 X 10-^ 6.40 x 10-^ 5.99 x 10"^ 
10 8.99 X 10^4 1.24 X lO'^ 1.08 x 10"^ 

3.3 Heat Model 



The full-order model is a plate with two heat sources and two points of measurements, and 
described by the heat equation as explained in [HI [IE]- A model of order n = 197 is obtained 
by spatial discretization. We choose a SISO subsystem corresponding to the first input and 
first output. Using IHA, we reduce the order to r = 2,4 and r = 6. As in the previous 
examples, we first tabulate, in Table |5} the results for solving the Step 2 of Algorithm IHA 
exactly and by the method of §2.2[ The conclusion is the same as before: The Loewner 
matrix approach yields T-Loo error norms and optimal dr values which are very close to the 
true optimal values of the underlying optimization problem in Step 2 of Algorithm IHA, indeed 
exact to the fifth digit for the r = 6 case. The decay of the singular values of SjL — M is 
shown in Figure [6} illustrating the choice of k. Only the r = Q case is presented; the other 
cases show the same pattern. 
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Figure 5: T/qo error as a function of the dr-term in MBT and IHA 



Table 5: Solution of the optimization problem in Step 2 
Exact Loewner 



dl 



\H — Ht 



IWc 



dl 



\H — Ht 



I lio 



k 



-2.0694 X 10-1 1.0710 x 10"^ 

2.0875 X 10-2 8.9082 x 10"^ 

-5.6250 X 10~^ 2.3578 x 10"^ 



-2.0700 X 10-1 

2.0813 X 10-2 

-5.6250 X 10-1 



1.0711 X 10-2 7 
8.9166 X 10-^ 9 
2.3578 X 10-^ 13 



Results for comparison with MBT and HNA are shown in Table |6} Once more, the proposed 
method consistently yields better 1-L^ performance than MBT. Even though for r = 2 the 
proposed method leads to smaller "Hoc error, for r = 4, 6, HNA yields slightly better re- 
sults. Hence, for this example as well, using interpolatory projections, we are able to beat 
MBT consistently and yield results comparable to or better than that of HNA. Indeed, this is 
satisfactory since, as stated in the introduction, for the large-scale settings we are interested 
in, implementing HNA will be a formidable task, if not impossible. 

As done in the previous examples, we illustrate, in Figure [7], the behavior of the absolute 
?^oo-error while optimizing over the (i^-term in both IHA and MBT. In this case, MBT almost 
gains nothing from the (i^-term modification, as the T^oo-error is reduced from 1.0897 x 10"'^ 
only to 1.0894 x 10~^, a marginal gain of 0.027%. On the other hand, IHA reduces the T-Lca- 
error from 1.26 x 10-'^ down to 5.46 x 10^^, a reduction factor of 56.89%. Once again, this 
reduction in the error is achieved even when the initial Hoo-error is bigger than that of BT, 
providing a clear illustration of the effectiveness of the (ir-term optimization in the proposed 
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Normalized singular values of s L — M for r=6 




Figure 6: The decay of the singular values of SjL — M for the Heat Model 

Table 6: Relative "Hoc error norms for the Heat Model 

r IHA MBT HNA 

~2 1.08 X 10^ 1.66 X 10-2 1.11 X 10"^ 

4 8.92 X 10-^ 1.68 x lO^^ 8.47 x 10"^ 

6 2.30 X 10-5 4.61 X 10"^ 2.07 x 10"^ 

method. 

To illustrate the importance of using the IRKA points in Step 2 of Algorithm IHA, we use 
arbitrarily chosen interpolation points (within the bounds of the mirror spectrum of the full 
order A) rather than the IRKA points. Then, we apply the same (i^-term modification as 
before. For r = 6, for example, the resulting relative T^oo-error is 1.72 x 10"^, two order 
of magnitudes higher than what we obtain using the IRKA points. This simple example 
illustrates the advantage of initializing Step 2 of Algorithm IHA with points computed by 
IRKA. We want to emphasize that these arbitrary interpolation points are indeed used to 
initialize IRKA. Hence, IRKA corrects these arbitrarily chosen points, producing interpolation 
points that are used to obtain an T^oo-error norm two orders of magnitude smaller. 



3.4 A Heat Transfer Problem in the Cooling of Steel Profiles 

This model describes a cooling process in a rolling mill which has been modeled as boundary 
control of a two dimensional heat equation. The full order model has n = 79841 with 7 
inputs and 6 outputs. We consider a SISO subsystem corresponding to the sixth input and 
the second output. For details regarding this model, see [8l|9]. For such a large order, imple- 
menting HNA is not possible. BT can be implemented iteratively using ADI-type methods; 
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Figure 7: T/qo error as a function of the dr-term in MBT and IHA 



however this requires state-of-the-art iterative Lyapunov solvers (two generahzed Lyapunov 
equations of order n = 79841 need to be solved) and is not the focus of this paper. Hence, we 
concentrate only on the performance of the proposed method and compare it with IRKA to 
show the improvement by the dr-terra optimzation. IRKA in Step 1 of Algorithm IHA is im- 
plemented in Matlab using direct sparse linear solves. Once again, we consider the solution 
of the optimization problem in Step 2 of Algorithm IHA using both approaches, i.e., directly 
solving the large-scale optimization problem versus using the method of §2.2 The results are 
shown in Table [7] and reveal the same pattern as before. Note that the direct solution of this 
scalar optimization problem requires several Hoo norm computation for a system of order 
79841 + r. This is computationally intractable, so the Tioo norms for the direct method are 
computed approximately by sampling along the imaginary axis at 500 points, logarithmically 
spaced points between 10~^ and 10. This is not an issue for the Loewner matrix approach 
since T^oo-norm computations are done on surrogate error systems of order k; which does not 
exceed 13 for this example. The decay of the singular values of SjL — M is shown in Figure 



Table 7: Solution of the optimization problem in Step 2 





Exact 




Loewner 




r 


d: 


\\H - ^rW-Hoc 


rf* \H-H:\n^ 


k 


2 


1.0189 X 10-2 


6.3715 X 10-1 


1.2685 X 10-2 6.3725 x 10-^ 


7 


4 


-1.0500 X 10"^ 


7.4620 X 10-2 


-5.3232 X 10-5 7.5483 x 10-^ 


12 


6 


-1.6859 X 10-^ 


5.4567 X 10-=^ 


-1.6849 X 10-5 5.4592 x 10-^ 


13 
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[8} indicating a natural choice for k. Only the r = 6 case is presented; the other cases show 
the same pattern. 



Normalized singular values of s L — M for r=6 




Figure 8: The decay of the singular values of SjL — M for the Cooling Steel Model 

We compare the performance of IHA with that of IRKA for r = 2,4,6, to illustrate the 
improvement offered by optimizing over the rf^-term. The results are listed in Table [8] As 
expected, IHA outperforms IRKA for every r- value. 

Table 8: Relative T-Loo error for the Cooling Steel Model, order = 79841 



IRKA 



IHA 



2 6.46 X 10-1 6.37 x lO^^ 
4 1.80 X 10-1 7.46 X lO^^ 
6 1.40 X 10-2 5.46 X 10"^ 



In Figure |9j we demonstrate how the absolute Tioo error changes over values of dr, for the 
order r = 6 approximation. In this case, the "Hoo-error is decreased by over a factor of two, 
from 2.58 x 10~^ down to 1.09 x 10~^. We also note that, for r = 6, computing the optimal 
d* term placed an additional interpolation point at 2.44 x 10"^, yielding exactly 2r + 1 = 13 
interpolation points in C_|. as suggested by Theorem 2.1 Recalling Theorem |2.1 in Figure 



10 we demonstrate how this additional interpolation condition therefore results in a tighter. 



more circular Nyquist plot, which is equivalent to the image of the error along the imaginary 
axis being nearly circular. 
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iiorrn of tlie error systein as d,. varies 




Figure 9: T-ioo error as a function of the dr-ieiTu 

4 Conclusions 

We have introduced an interpolation-based model reduction technique to construct high- 
fidelity T/oo approximations for large-scale linear dynamical systems. For a given order r, r 
Hermite interpolation points are produced that induce (locally) optimal I-L2 approximation. 
The dr (feed-forward) term is then adjusted in a such way that interpolation at these initial 
2r points is retained while an additional interpolation point is added that minimizes the "Hoo- 
error norm. By employing a data-driven Loewner approach, this last step may be performed 
at negligible cost; no large-scale "Hoo-norm computations are ever needed. The dominant 
cost of the method lies with the solution of sparse linear systems. Four numerical examples 
show that the proposed method produces high fidelity "Hoo reduced-order models that are 
better than those obtained by balanced truncation; and as good as (and sometimes better 
than) those obtained by optimal Hankel norm approximation; in all cases, at much lower 
computational cost. 
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